# 
# --- Analysis file for "Free Trade and Forms of Democratization" (The Journal of Politics).
# --- This file calculates and exports virtually all tables and figures in the main manuscript and appendices.*
# --- It requires the trains2.csv and trains2_countrylevel.csv files provided in the replication package.
# --- Variable descriptions are provided under the "DATA LOADING" header.
#
# --- *The analyses in appendix E are done in the Stata .do file included in the replication package.
#
# --- Contact: noah.zucker@columbia.edu
#
# --- Last edited April 27, 2020.
#

rm(list = ls())

pacman::p_load(tidyverse, # for tidyverse operations
               broom, # for tidying looped regression results
               ggpubr, # for organizing ggplot graphics
               lfe, # for regression analyses
               stargazer, # for summary statistics and regression tables
               xtable) # for transition list table

#################### #
#### DATA LOADING ####
#################### #

# Note that lagged variables include suffixes denoting lag length in years (_l1, _l3, _l5, _l10)
#
# ccode: IMF country code
# country: country name
# year: year
# hs2: two-digit Harmonized System code
# tariff_smp_mfn: MFN tariff rate, simple mean at HS2 level
# polity2: Polity2 score
# mass: whether a country experienced a mass-led transition in a given year (overlaps with elite)
# elite: whether a country experienced an elite-led transition in a given year (overlaps with mass)
# mass2: whether a country experienced a mass-led transition in a given year (no overlaps with elite)
# elite2: whether a country experienced an elite-led transition in a given year (no overlaps with mass)
# gdp_per_capita_constant_2010_us: GDP per capita in constant 2010 US dollars
# imports_of_goods_and_services_percent_of_gdp: imports of goods/services as % of GDP
# population_total: national population
# use_of_imf_credit_dod_current_us: use of IMF credit, amount in current US dollars
# WTO: indicator for WTO membership (= 1)
# idealpoint: UNGA voting ideal point of country (Bailey, Strehnev, and Voeten; JCR 2017)
# idealpoint_US: UNGA voting ideal point of US (Bailey, Strehnev, and Voeten; JCR 2017)
# idealpoint_diff: difference in above UNGA voting ideal points (Bailey, Strehnev, and Voeten; JCR 2017)
# oda_total_gross_disbursements_dac: total gross ODA disbursements to country from DAC members
# rca: revealed comparative advantage for given product group

trains2 <- read_csv("trains2.csv", col_types = cols(.default = "d",
                                                    country = "c",
                                                    transition_polity_year = "c"))
trains2c <- read_csv("trains2_countrylevel.csv", col_types = cols(.default = "d",
                                                                  country = "c",
                                                                  transition_polity_year = "c"))

############################ #
#### DESCRIPTIVE ANALYSES ####
############################ #

# > Before/after density plots of tariff rates (figures A1-A2) ####

.mass_countries2 <- unique(trains2$ccode[trains2$mass == 1])[-1]
.elite_countries2 <- unique(trains2$ccode[trains2$elite == 1])[-1]

trains3 <- trains2 %>% filter(!is.na(mass) | !is.na(elite))

for (i in 1:(nrow(trains3) - 1)){
  trains3$mass[i] <- ifelse(trains3$ccode[i+1] == trains3$ccode[i] &
                              trains3$year[i+1] > trains3$year[i] &
                              trains3$mass[i+1] == 1 & trains3$mass[i] == 0,
                            -1, trains3$mass[i])
  trains3$elite[i] <- ifelse(trains3$ccode[i+1] == trains3$ccode[i] &
                               trains3$year[i+1] > trains3$year[i] &
                               trains3$elite[i+1] == 1 & trains3$elite[i] == 0,
                            -1, trains3$elite[i])
}

trains3 %>%
  filter(mass == -1) %>%
  filter(!is.na(ccode)) %>%
  select(ccode, hs2,tariff_smp_mfn) %>%
  mutate(type_group = "Mass") -> temp_pre.mass.transition

trains3 %>%
  filter(elite == -1) %>%
  filter(!is.na(ccode)) %>%
  select(ccode, hs2,tariff_smp_mfn) %>%
  mutate(type_group = "Elite") -> temp_pre.elite.transition

temp_pre.transition <- bind_rows(temp_pre.elite.transition, temp_pre.mass.transition)

trains2 %>%
  filter(!is.na(ccode)) %>%
  mutate(regime_type = case_when(polity2 >= 6 ~ "Democracy",
                                 polity2 < 6  ~ "Non-democracy",
                                 TRUE ~ "other")) %>%
  filter(regime_type != "other") %>%
  ggplot() +
  geom_density(aes(x = log1p(tariff_smp_mfn), 
                   group = regime_type, 
                   fill = regime_type), 
               alpha = .4) +
  theme_minimal() +
  ylim(0,1) +
  xlim(0,7) +
  xlab("HS2-level tariff rates (ln)") +
  scale_fill_viridis_d("Regime Type (Polity2)",
                       option = "C") +
  theme(axis.line = element_blank(),
        axis.ticks.y = element_blank(),        
        axis.title.x = element_text(size=9),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.position = "bottom",
        plot.title = element_text(size=12, face = "bold"),
        legend.title = element_text(size=9, face = "bold")) -> dem_autoc_plot

temp_pre.transition %>%
  ggplot() +
  geom_density(aes(x = log1p(tariff_smp_mfn), 
                   group = type_group, 
                   fill = type_group), 
               alpha = .4) +
  theme_minimal() +
  ylim(0,1) +
  xlim(0,7) +
  xlab("HS2-level tariff rate in year before transition (ln)") +
  scale_fill_viridis_d("Transition Type") +
  theme(axis.line = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_text(size=9),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.position = "bottom",
        plot.title = element_text(size=12, face = "bold"),
        legend.title = element_text(size=9, face = "bold")) -> trains2_pre.transition_plot

trains2 %>%
  mutate(type_group = NA) %>%
  mutate(type_group = replace(type_group, elite_l1 == 1, "Elite"),
         type_group = replace(type_group, mass_l1 == 1, "Mass")) %>%
  filter(!is.na(type_group)) %>%
  ggplot() +
  geom_density(aes(x = log1p(tariff_smp_mfn), 
                   group = type_group, 
                   fill = type_group), 
               alpha = .4) +
  theme_minimal() +
  ylim(0,1) +
  xlim(0,7) +
  xlab("HS2-level tariff rates after transition (ln)") +
  scale_fill_viridis_d("Transition Type") +
  theme(axis.line = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_text(size=9),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.position = "bottom",
        plot.title = element_text(size=12, face = "bold"),
        legend.title = element_text(size=9, face = "bold")) -> trains2_post.transition_plot

ggarrange(trains2_pre.transition_plot, trains2_post.transition_plot,
          ncol = 2, 
          common.legend = TRUE, legend = "bottom") -> density.plots2

ggsave(filename = "densityplots2.pdf",
       plot = density.plots2,
       width = 8, height = 4)

ggsave(filename = "dem_autoc_plot.pdf",
       plot = dem_autoc_plot,
       width = 4, height = 4)

# > Summary statistics (table A1) ####

trains2 %>%
  as.data.frame() %>%
  select(tariff_smp_mfn,
         year,
         mass2,
         elite2,
         mass,
         elite,
         rca,
         gdp_per_capita_constant_2010_us,
         WTO,
         imports_of_goods_and_services_percent_of_gdp,
         population_total,
         oda_total_gross_disbursements_dac,
         use_of_imf_credit_dod_current_us,
         idealpoint_diff) %>%
  stargazer(covariate.labels = c("AVE MFN tariff rates, HS2 level",
                                 "Year",
                                 "Elite-led transition",
                                 "Mass-led transition",
                                 "Elite-led transition (w/ overlap)",
                                 "Mass-led transition (w/ overlap)",
                                 "Product group RCA",
                                 "GDP/capita",
                                 "Full WTO membership",
                                 "Imports (\\% GDP)",
                                 "Total population",
                                 "Gross ODA received",
                                 "Use of IMF credit",
                                 "Distance from U.S. and UNGA"),
            float = F,
            no.space = T,
            omit.summary.stat = c("p25", "p75"),
            digits = 0,
            out = "summary_trains2.tex")

# > Transition list (table A3) ####

trains2 %>%
 filter(mass == 1) %>%
  distinct(country.x, mass, transition_polity_year) %>%
  rename(Country = country.x) %>%
  rename("Type" = mass) %>%
  mutate(Type = "Mass-led") %>%
  rename("Year" = transition_polity_year) %>%
  arrange(Country) -> mass_list

trains2 %>%
  filter(elite == 1) %>%
  distinct(country.x, elite, transition_polity_year) %>%
  rename(Country = country.x) %>%
  rename("Type" = elite) %>%
  mutate(Type = "Elite-led") %>%
  rename("Year" = transition_polity_year) %>%
  arrange(Country) -> elite_list

full_list <- bind_rows(mass_list, elite_list)

# Note: need to check and manually adjust countries with multiple transitions (e.g., Zambia)
print(xtable(full_list, type = "latex",
                     digits=c(0,0,0,0)), 
      include.rownames = F, 
      floating = F,
      tabular.environment = "longtable")

########################################################### #
#### MAIN REGRESSION MODELS (table 1 in main manuscript) ####
########################################################### #

# > With RCA ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_e1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_m1

# > All covariates ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1)
     | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m1

# > Exporting ####

stargazer(felm_rca_e1, felm_rca_m1, felm_cov_e1, felm_cov_m1,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_main.tex")

####################################### #
#### SUPPLEMENTARY REGRESSION MODELS ####
####################################### #

# > By-HS2 models (figure 1) ####

felm_by_hs2 <- vector("list", n_distinct(trains2$hs2))
hs2_list <- unique(trains2$hs2)
felm_by_hs2_results <- tibble(HS2 = sort(hs2_list),
                              Coefficient = NA,
                              SE = NA,
                              p = NA)

for (i in 1:n_distinct(trains2$hs2)){
  felm_by_hs2[[i]] <- felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
                             log(gdp_per_capita_constant_2010_us_l1) +
                             as.numeric(WTO_l1 == 1) +
                             imports_of_goods_and_services_percent_of_gdp_l1 +
                             log(population_total_l1) | ccode + year | 0 | ccode + year, 
                           data = trains2[trains2$hs2 == hs2_list[[i]], ],
                           exactDOF = T)
  
  felm_by_hs2_results$Coefficient[i] <- tidy(felm_by_hs2[[i]])$estimate[1]
  felm_by_hs2_results$SE[i] <- tidy(felm_by_hs2[[i]])$std.error[1]
  felm_by_hs2_results$p[i] <- tidy(felm_by_hs2[[i]])$p.value[1]
}

felm_by_hs2_results %>%
  ggplot() +
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = HS2, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), 
                 col = "gray80") +
  geom_point(aes(x = HS2, y = Coefficient, col = ifelse(p < .05, "black", "gray80"))) +
  scale_color_manual(values = c("black", "gray80")) +
  lims(y = c(-1, 1)) +
  scale_x_continuous(breaks = c(1, seq(10,90,10), 97), labels = c(1, seq(10,90,10), 97)) +
  theme_minimal() +
  ggtitle("(a) Elite-led transitions") +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank()) -> felm_by_hs2_results_plot

felm_by_hs2_mass <- vector("list", n_distinct(trains2$hs2))
hs2_list <- unique(trains2$hs2)
felm_by_hs2_mass_results <- tibble(HS2 = sort(hs2_list),
                                   Coefficient = NA,
                                   SE = NA,
                                   p = NA)

for (i in 1:n_distinct(trains2$hs2)){
  felm_by_hs2_mass[[i]] <- felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
                                  log(gdp_per_capita_constant_2010_us_l1) +
                                  as.numeric(WTO_l1 == 1) +
                                  imports_of_goods_and_services_percent_of_gdp_l1 +
                                  log(population_total_l1) | ccode + year | 0 | ccode + year, 
                                data = trains2[trains2$hs2 == hs2_list[[i]], ],
                                exactDOF = T)
  
  felm_by_hs2_mass_results$Coefficient[i] <- tidy(felm_by_hs2_mass[[i]])$estimate[1]
  felm_by_hs2_mass_results$SE[i] <- tidy(felm_by_hs2_mass[[i]])$std.error[1]
  felm_by_hs2_mass_results$p[i] <- tidy(felm_by_hs2_mass[[i]])$p.value[1]
}

felm_by_hs2_mass_results %>%
  ggplot() +
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = HS2, ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE), 
                 col = "gray80") +
  geom_point(aes(x = HS2, y = Coefficient, col = ifelse(p < .05, "black", "gray80"))) +
  scale_color_manual(values = c("black", "gray80")) +
  lims(y = c(-2, 1)) +
  scale_x_continuous(breaks = c(1, seq(10,90,10), 97), labels = c(1, seq(10,90,10), 97)) +
  ggtitle("(b) Mass-led transitions") +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank()) -> felm_by_hs2_mass_results_plot

ggarrange(felm_by_hs2_results_plot, felm_by_hs2_mass_results_plot,
                  nrow = 2) -> felm_by_hs2_results_plot_combined

ggsave(felm_by_hs2_results_plot_combined,
       device = "eps",
       dpi = 300,
       filename = "felm_by_hs2_results_plot_combined.eps",
       width = 6.5,
       height = 4.5,
       units = "in")

# > By-ccode models (figure D1) ####

ccode_list <- unique(trains2$ccode)
felm_by_ccode <- vector("list", length(ccode_list))
felm_by_ccode_results <- tibble(ccode = sort(ccode_list),
                                Coefficient = NA,
                                SE = NA,
                                p = NA,
                                df = NA)

for (i in 1:n_distinct(trains2$ccode)){
  .trains2_temp <- trains2[trains2$ccode != ccode_list[[i]],]
  
  felm_by_ccode[[i]] <- felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
                               log(gdp_per_capita_constant_2010_us_l1) +
                               as.numeric(WTO_l1 == 1) +
                               imports_of_goods_and_services_percent_of_gdp_l1 +
                               log(population_total_l1) | ccode + year | 0 | ccode + year, 
                             data = .trains2_temp,
                             exactDOF = T)
  
  felm_by_ccode_results$Coefficient[i] <- tidy(felm_by_ccode[[i]])$estimate[1]
  felm_by_ccode_results$SE[i] <- tidy(felm_by_ccode[[i]])$std.error[1]
  felm_by_ccode_results$p[i] <- tidy(felm_by_ccode[[i]])$p.value[1]
  felm_by_ccode_results$df[i] <- summary(felm_by_ccode[[i]])$df[1]
}

felm_by_ccode_results %>%
  ggplot() +
  geom_point(aes(x = 1:nrow(felm_by_ccode_results), y = Coefficient)) +
  geom_linerange(aes(x = 1:nrow(felm_by_ccode_results), 
                     ymin = Coefficient - 1.96*SE, ymax = Coefficient + 1.96*SE),
                 alpha = .3) +
  geom_hline(yintercept = 0, lty = 2) +
  lims(y = c(-.5, .5)) +
  scale_x_continuous(breaks = c(1, 42, 84, 128, 168), labels = c(111, 361, 578, 742, 968)) +
  theme_minimal() +
  labs(x = "IMF Country Code") +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank()) -> felm_by_ccode_results_plot

ggsave(felm_by_ccode_results_plot,
       filename = "iterative_sampling2.pdf",
       width = 6.5,
       height = 2,
       units = "in")

# > Conditioning on pre-transition rates (table D1) ####

temp_pre.mass.transition <- temp_pre.mass.transition %>%
  select(-type_group) %>%
  rename(tariff_smp_mfn_pre.mass = tariff_smp_mfn)

temp_pre.elite.transition <- temp_pre.elite.transition %>%
  select(-type_group) %>%
  rename(tariff_smp_mfn_pre.elite = tariff_smp_mfn)

trains4 <- trains2 %>%
  left_join(temp_pre.mass.transition, by = c("ccode", "hs2")) %>%
  left_join(temp_pre.elite.transition, by = c("ccode", "hs2"))

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + log1p(tariff_smp_mfn_pre.elite) | ccode + year | 0 | ccode + year, 
     data = trains4,
     exactDOF = T) -> felm_rca_e1sp1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + log1p(tariff_smp_mfn_pre.mass) | ccode + year | 0 | ccode + year, 
     data = trains4,
     exactDOF = T) -> felm_rca_m1sp1

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) + 
       log1p(tariff_smp_mfn_pre.elite) | ccode + year | 0 | ccode + year, 
     data = trains4,
     exactDOF = T) -> felm_cov_e1sp2

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) + 
       log1p(tariff_smp_mfn_pre.mass)
     | ccode + year | 0 | ccode + year, 
     data = trains4,
     exactDOF = T) -> felm_cov_m1sp2

stargazer(felm_rca_e1sp1, felm_rca_m1sp1, felm_cov_e1sp2, felm_cov_m1sp2,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)",
                               "Pre-elite transition tariff rates (ln)",
                               "Pre-mass transition tariff rates (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_startingpoints.tex")

# > Three-year lags (table D2) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l3 + rca_l3 + 
       log(gdp_per_capita_constant_2010_us_l3) +
       as.numeric(WTO_l3 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l3 +
       log(population_total_l3) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e3

felm(log1p(tariff_smp_mfn) ~ mass2_l3 + rca_l3 + 
       log(gdp_per_capita_constant_2010_us_l3) +
       as.numeric(WTO_l3 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l3 +
       log(population_total_l3) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m3

stargazer(felm_cov_e3, felm_cov_m3,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_lag3.tex")

# > Five-year lags (table D3) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l5 + rca_l5 + 
       log(gdp_per_capita_constant_2010_us_l5) +
       as.numeric(WTO_l5 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l5 +
       log(population_total_l5) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e5

felm(log1p(tariff_smp_mfn) ~ mass2_l5 + rca_l5 + 
       log(gdp_per_capita_constant_2010_us_l5) +
       as.numeric(WTO_l5 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l5 +
       log(population_total_l5) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m5

stargazer(felm_cov_e5, felm_cov_m5,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_lag5.tex")

# > Ten-year lags (table D4) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l10 + rca_l10 + 
       log(gdp_per_capita_constant_2010_us_l10) +
       as.numeric(WTO_l10 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l10 +
       log(population_total_l10) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e10

felm(log1p(tariff_smp_mfn) ~ mass2_l10 + rca_l10 + 
       log(gdp_per_capita_constant_2010_us_l10) +
       as.numeric(WTO_l10 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l10 +
       log(population_total_l10) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m10

stargazer(felm_cov_e10, felm_cov_m10,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_lag10.tex")

# > Elite, mass, polity togther; elite and mass together (table D5) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + mass2_l1 + polity2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_emp

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_emp2

stargazer(felm_cov_emp2,
          felm_cov_emp,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Polity2",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_emp.tex")

# Wald test (caption of table D5)
lfe::waldtest(felm_cov_emp2, R = ~ mass2_l1 - elite2_l1)

# > Overlapping transition indicators (table D6) ####

felm(log1p(tariff_smp_mfn) ~ elite_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_e1o

felm(log1p(tariff_smp_mfn) ~ mass_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_m1o

felm(log1p(tariff_smp_mfn) ~ elite_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e1o

felm(log1p(tariff_smp_mfn) ~ mass_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1)
     | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m1o

stargazer(felm_rca_e1o, felm_rca_m1o, felm_cov_e1o, felm_cov_m1o,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition (w/ overlap)",
                               "Mass-led transition (w/ overlap)",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_nonas.tex")

# > Clustering only by country (table D7) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode, 
     data = trains2,
     exactDOF = T) -> felm_rca_e1c1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode, 
     data = trains2,
     exactDOF = T) -> felm_rca_m1c1

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode, 
     data = trains2,
     exactDOF = T) -> felm_cov_e1c1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1)
     | hs2 + ccode + year | 0 | ccode, 
     data = trains2,
     exactDOF = T) -> felm_cov_m1c1

stargazer(felm_rca_e1c1, felm_rca_m1c1, felm_cov_e1c1, felm_cov_m1c1,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_onewayclustering.tex")

# > Clustering by country, industry, and year (table D8) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 | hs2 + ccode + year | 0 | hs2 + ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_e1c3

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 | hs2 + ccode + year | 0 | hs2 + ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_m1c3

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | hs2 + ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e1c3

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1)
     | hs2 + ccode + year | 0 | hs2 + ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m1c3

stargazer(felm_rca_e1c3, felm_rca_m1c3, felm_cov_e1c3, felm_cov_m1c3,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_threewayclustering.tex")

# > Country and year fixed effects only (table D9) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 | ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_e1fe2

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 | ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_rca_m1fe2

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_e1fe2

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1)
     | ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_cov_m1fe2

stargazer(felm_rca_e1fe2, felm_rca_m1fe2, felm_cov_e1fe2, felm_cov_m1fe2,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_twowayfixedeffects.tex")

# > Post-1995 sample (table D10) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2[trains2$year >= 1995,],
     exactDOF = T) -> felm_post95_e1

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2[trains2$year >= 1995,],
     exactDOF = T) -> felm_post95_m1

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2[trains2$year >= 1995,],
     exactDOF = T) -> felm_post95_e2

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2[trains2$year >= 1995,],
     exactDOF = T) -> felm_post95_m2

stargazer(felm_post95_e1, felm_post95_m1, felm_post95_e2, felm_post95_m2,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_post95.tex")

# > Foreign pressure covariates (table D11) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) +
       log1p(oda_total_gross_disbursements_dac_l1) +
       log1p(use_of_imf_credit_dod_current_us_l1) +
       idealpoint_diff_l1
     | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_foreign_e2

felm(log1p(tariff_smp_mfn) ~ mass2_l1 + rca_l1 + 
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) +
       log1p(oda_total_gross_disbursements_dac_l1) +
       log1p(use_of_imf_credit_dod_current_us_l1) +
       idealpoint_diff_l1
     | hs2 + ccode + year | 0 | ccode + year, 
     data = trains2,
     exactDOF = T) -> felm_foreign_m2

stargazer(felm_foreign_e2, felm_foreign_m2,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)",
                               "Gross ODA received (ln)",
                               "Use of IMF credit (ln)",
                               "Distance from U.S. at UNGA"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_foreign.tex")

# > Lagged DV models (table D12) ####

trains2 %>%
  group_by(ccode, hs2) %>%
  mutate(tariff_smp_mfn_l1 = dplyr::lag(tariff_smp_mfn)) %>%
  ungroup() %>%
  felm(log1p(tariff_smp_mfn) ~ elite2_l1 + log1p(tariff_smp_mfn_l1) | hs2 + ccode + year | 0 | ccode + year, 
       data = .,
       exactDOF = T) -> felm_laggeddv_e1

trains2 %>%
  group_by(ccode, hs2) %>%
  mutate(tariff_smp_mfn_l1 = dplyr::lag(tariff_smp_mfn)) %>%
  ungroup() %>%
  felm(log1p(tariff_smp_mfn) ~ mass2_l1 + log1p(tariff_smp_mfn_l1) | hs2 + ccode + year | 0 | ccode + year, 
       data = .,
       exactDOF = T) -> felm_laggeddv_m1

trains2 %>%
  group_by(ccode, hs2) %>%
  mutate(tariff_smp_mfn_l1 = dplyr::lag(tariff_smp_mfn)) %>%
  ungroup() %>%
  felm(log1p(tariff_smp_mfn) ~ elite2_l1 + log1p(tariff_smp_mfn_l1) + rca_l1 + 
         log(gdp_per_capita_constant_2010_us_l1) +
         as.numeric(WTO_l1 == 1) +
         imports_of_goods_and_services_percent_of_gdp_l1 +
         log(population_total_l1) | hs2 + ccode + year | 0 | ccode + year, 
       data = .,
       exactDOF = T) -> felm_laggeddv_e2

trains2 %>%
  group_by(ccode, hs2) %>%
  mutate(tariff_smp_mfn_l1 = dplyr::lag(tariff_smp_mfn)) %>%
  ungroup() %>%
  felm(log1p(tariff_smp_mfn) ~ mass2_l1 + log1p(tariff_smp_mfn_l1) + rca_l1 + 
         log(gdp_per_capita_constant_2010_us_l1) +
         as.numeric(WTO_l1 == 1) +
         imports_of_goods_and_services_percent_of_gdp_l1 +
         log(population_total_l1)
       | hs2 + ccode + year | 0 | ccode + year, 
       data = .,
       exactDOF = T) -> felm_laggeddv_m2

stargazer(felm_laggeddv_e1, felm_laggeddv_m1, felm_laggeddv_e2, felm_laggeddv_m2,
          dep.var.labels = "AVE MFN tariff rates, HS2 level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "Lagged DV (tariff rates)",
                               "Product group RCA",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_laggeddv.tex")

# > Country-level regressions (table D13) ####

felm(log1p(tariff_smp_mfn) ~ elite2_l1 | ccode + year | 0 | ccode + year, 
     data = trains2c,
     exactDOF = T) -> felm_bv_e1cl

felm(log1p(tariff_smp_mfn) ~ mass2_l1 | ccode + year | 0 | ccode + year, 
     data = trains2c,
     exactDOF = T) -> felm_bv_m1cl

felm(log1p(tariff_smp_mfn) ~ elite2_l1 +
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | ccode + year | 0 | ccode + year, 
     data = trains2c,
     exactDOF = T) -> felm_cov_e1cl

felm(log1p(tariff_smp_mfn) ~ mass2_l1 +
       log(gdp_per_capita_constant_2010_us_l1) +
       as.numeric(WTO_l1 == 1) +
       imports_of_goods_and_services_percent_of_gdp_l1 +
       log(population_total_l1) | ccode + year | 0 | ccode + year, 
     data = trains2c,
     exactDOF = T) -> felm_cov_m1cl

stargazer(felm_bv_e1cl, felm_bv_m1cl, 
          felm_cov_e1cl, felm_cov_m1cl,
          dep.var.labels = "AVE MFN tariff rates, country-year level (ln)",
          covariate.labels = c("Elite-led transition",
                               "Mass-led transition",
                               "GDP/capita (ln)",
                               "Full WTO membership",
                               "Imports (\\% GDP)",
                               "Total population (ln)"),
          float = F,
          omit.stat = c("f", "ser", "rsq"),
          no.space = T,
          out = "felm_cl.tex")

